module mRoot

        use mPrecision
        use mObject

        implicit none

        interface new
                module procedure new_root
        end interface new

        interface getParam
                module procedure getParam_root
        end interface getParam

        interface getParams
                module procedure getParams_2_root
        end interface getParams

        interface getLPS
                module procedure getLPS_root
        end interface getLPS

        interface getChi
                module procedure getChi_root
        end interface getChi

        interface getSmith
                module procedure getSmith_root
        end interface getSmith

        private
        public :: new, findRoot, getParam, getParams, getLPS, getChi, getSmith, getOutput, setOutput

contains

        subroutine new_root(root, ll, ul, tolerance, params, lps, chi, smith)

                use mMisc

                type(tRoot), intent(out) :: root
                real(D), intent(in) :: ll, ul
                real(D), intent(in), optional :: tolerance, params(:)
                type(tLPS), intent(in), optional :: lps
                type(tChi), intent(in), optional :: chi
                type(tSmith), intent(in), optional :: smith

                real(D), parameter :: DEFAULT_TOLERANCE = epsilon(1.0_D)

                root%ll = ll
                root%ul = ul

                if( present(tolerance) ) then
                        root%tolerance = tolerance
                else
                        root%tolerance = DEFAULT_TOLERANCE
                end if

                if( present(params) ) then
                        call assert( size(params) <= ROOT_N_PARAMS, 'Too many params to root' )
                        root%params( 1:size(params) ) = params(:)
                end if

                if( present(lps) ) root%lps = lps
                if( present(chi) ) root%chi = chi
                if( present(smith) ) root%smith = smith
        end subroutine new_root

        real(D) function getParam_root(root)

                type(tRoot), intent(in) :: root

                getParam_root = root%params(1)

        end function getParam_root

        subroutine getParams_2_root(root, p1, p2)

                type(tRoot), intent(in) :: root
                real(D), intent(out) :: p1, p2
                
                p1 = root%params(1)
                p2 = root%params(2)

        end subroutine getParams_2_root

        type(tLPS) function getLPS_root(root)
                type(tRoot), intent(in) :: root

                getLPS_root = root%lps
        end function getLPS_root

        type(tChi) function getChi_root(root)
                type(tRoot), intent(in) :: root

                getChi_root = root%chi
        end function getChi_root

        type(tSmith) function getSmith_root(root)
                type(tRoot), intent(in) :: root

                getSmith_root = root%smith
        end function getSmith_root

        real(D) function getOutput(root)

                type(tRoot), intent(in) :: root

                getOutput = root%output

        end function getOutput

        subroutine setOutput(root, val)

                type(tRoot), intent(inout) :: root
                real(D), intent(in) :: val

                root%output = val
        end subroutine setOutput

        !********************************
        ! Brent's Method for root finding
        !********************************
!        subroutine findRoot(root, func, ans, out)
        logical function findRoot(root, func, ans, out)

                type(tRoot), intent(inout) :: root
                interface
                        function func(root, x)
                                use mPrecision
                                use mObject
                                real(D) :: func
                                type(tRoot), intent(inout) :: root
                                real(D), intent(in) :: x
                        end function func
                end interface
                real(D), intent(out) :: ans
                real(D), intent(out), optional :: out

                integer, parameter :: ITMAX = 100
                real(D), parameter :: EPS = epsilon(1.0_D)

                integer :: iter
                real(D) :: a, b, c, dHere, e, fa, fb, fc, p, q, r, s, tol1, xm

                a = root%ll
                b = root%ul

                fa = func(root, a)
                fb = func(root, b)

                if( (fa .gt. 0.0_D .and. fb .gt. 0.0_D) &
                    .or. (fa .lt. 0.0_D .and. fb .lt. 0.0_D) ) then
                        !write(0, *) 'root must be bracketed for Brent. You have'
                        !write(0, *) fa, '<= root <=', fb
                        findRoot = .false.
                        return
                end if

                c = b
                fc = fb
                do iter = 1, ITMAX
                        if( (fb .gt. 0.0_D .and. fc .gt. 0.0_D) &
                            .or. (fb .lt. 0.0_D .and. fc .lt. 0.0_D) ) then
                                c = a
                                fc = fa
                                dHere = b - a
                                e = dHere
                        end if
                        if( abs(fc) .lt. abs(fb) ) then
                                a = b
                                b = c
                                c = a
                                fa = fb
                                fb = fc
                                fc = fa
                        end if
                        tol1 = 2.0_D * EPS * abs(b) + 0.5_D * root%tolerance
                        xm = 0.5_D * (c - b)
                        if( abs(xm) .le. tol1 .or. fb .eq. 0.0_D ) then
                                ans = b
                                findRoot = .true.
                                return
                        endif
                        if( abs(e) .ge. tol1 .and. abs(fa) .gt. abs(fb) ) then
                                s = fb / fa
                                if(a .eq. c) then
                                        p = 2.0_D * xm * s
                                        q = 1.0_D - s
                                else
                                        q = fa / fc
                                        r = fb / fc
                                        p = s * ( 2.0_D * xm * q * (q - r) &
                                                  - (b - a) * (r - 1.0_D) )
                                        q = (q - 1.0_D) * (r - 1.0_D) &
                                                        * (s - 1.0_D)
                                end if
                                if(p .gt. 0.0_D) q = -q
                                p = abs(p)
                                if(2.0_D * p &
                                    .lt. min( 3.0_D * xm * q - abs(tol1 * q), &
                                              abs(e * q) )) then
                                        e = dHere
                                        dHere = p / q
                                else
                                        dHere = xm
                                        e = dHere
                                end if
                        else
                                dHere = xm
                                e = dHere
                        end if
                        a = b
                        fa = fb
                        if( abs(dHere) .gt. tol1 ) then
                                b = b + dHere
                        else
                                b = b + sign(tol1, xm)
                        end if

                        fb = func(root, b)
                        if( present(out) ) out = root%output
                end do

                write(0, *) 'Brent exceeding maximum iterations'

                ans = b

                findRoot = .true.
        end function findRoot
!C  (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3.

end module mRoot
